#### clear working enviornment ####
rm(list = ls())

#### load required packages ####
library(tidyverse)
library(Hmisc)
library(stargazer)
library(marginaleffects)

#### load the replication data ####
df <- readRDS("data/dem-transitions-replication-data.rds") %>%
  glimpse()

#### calculate country-years, mid-years, initiation years, and baseline probabilities of outcomes ####
# set up a matrix to create Table 1 in the main text 
tab2 <- matrix(NA, 2, 3)

rownames(tab2) <- c("Pr(MID Participation)", "Pr(Fatal Participation)")

colnames(tab2) <- c("Full Sample", "Democracies", "Autocracies")

# full sample
tab2[1,1] <- round(mean(df$mid), 2)
tab2[2,1] <- round(mean(df$fatal), 2)

# democracies
dem <- df %>%
  filter(v2x_regime >= 2)
tab2[1,2] <- round(mean(dem$mid), 2)
tab2[2,2] <- round(mean(dem$fatal), 2)

# autocracies
aut <- df %>%
  filter(v2x_regime < 2)
tab2[1,3] <- round(mean(aut$mid), 2)
tab2[2,3] <- round(mean(aut$fatal), 2)

# code to make table in latex
latex(tab2, file = "", title = "Baseline Probabilities of MID Participation and Fatal Participation (1900-2014)", booktabs = TRUE, numeric.dollar = TRUE)

# chi-squared differences between democracy and autocracy for conflict measures
test_1 <- df %>%
  mutate(dem = ifelse(v2x_regime >= 2, 1, 0)) %>%
  select(dem, mid) %>%
  filter(dem >= 0) %>% 
  glimpse()

chisq.test(test_1$dem, test_1$mid)

test_2 <- df %>%
  mutate(dem = ifelse(v2x_regime >= 2, 1, 0)) %>%
  select(dem, fatal) %>%
  filter(dem >= 0) %>% 
  glimpse()

chisq.test(test_2$dem, test_2$fatal)

#### calculate benchmarks to compare magnitude of effect size ####

# estimate benchmark correlates of conflict onset using covariates
out1 <- glm(mid ~ ln_cinc + ln_pecpp + as.factor(terr_dispute) + mid_peace_yrs + mid_peace_yrs_2 + mid_peace_yrs_3, data = df, family = binomial(link = "logit"))
summary(out1)

out2 <- glm(fatal ~ ln_cinc + ln_pecpp + as.factor(terr_dispute) + fatal_peace_yrs + fatal_peace_yrs_2 + fatal_peace_yrs_3, data = df, family = binomial(link = "logit"))
summary(out2)

# create Table E.1 for the appendix
stargazer(out1, out2, style = "apsr", dep.var.labels = c("MID Participation", "MID Inititation"),
          covariate.labels = c("ln(CINC)", "ln(PEC per capita)", "Territorial Dispute", "All MIDs peace years", "All MIDs peace years$^2$", "All MIDs peace years$^3$", "Fatal MIDs peace years", "Fatal MIDs peace years$^2$", "Fatal MIDs peace years$^3$"))

# plot change in conflict across the interquartile range of capabilities
gg_1 <- avg_comparisons(out1, variables = list("ln_cinc" = "iqr")) %>%
  mutate(outcome = "All MIDs")
gg_2 <- avg_comparisons(out2, variables = list("ln_cinc" = "iqr")) %>%
  mutate(outcome = "Fatal MIDs")

gg <- rbind(gg_1, gg_2)

# create Figure E.1 for the appendix
ggplot(data = gg, aes(x = outcome, y = -estimate, ymin = -conf.low, ymax = -conf.high)) +
  geom_pointrange(linewidth = 1) +
  geom_hline(yintercept = 0, linetype = "dashed", linewidth = 1) +
  annotate("text", x = 1.1, y = -.121, label = ".12", size = 3.75) +
  annotate("text", x = 2.1, y = -.043, label = ".04", size = 3.75) +
  theme_minimal() +
  ylim(-.2, .1) +
  labs(x = "Outcome",
       y = "Estimated Change in Probability of Conflict")
ggsave("figures/meaningful-effect-comparison.pdf", width = 6.5, height = 4.5, units = "in", bg = "white")

